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ABSTRACT 

Supermassive black holes are believed to be ubiquitous at the centers of galaxies. Measuring their 
masses is extremely challenging yet essential for understanding their role in the formation and evolution 
of cosmic structure. We present a direct measurement of the mass of a black hole in an active galactic 
nucleus (Arp 151) based on the motion of the gas responsible for the broad emission lines. By analyzing 
and modeling spectroscopic and photometric time series, we find that the gas is well described by a 
disk or torus with an average radius of 3.99± 1.25 light days and an opening angle of 68.91^2 degrees, 
viewed at an inclination angle of 67.8 ±7.8 degrees (that is, closer to face-on than edge-on). The black 
hole mass is inferred to be io 6 - 51±0 - 28 M Q . The method is fully general and can be used to determine 
the masses of black holes at arbitrary distances, enabling studies of their evolution over cosmic time. 
Subject headings: galaxies: active — methods: data analysis — methods: statistical 



1. INTRODUCTION 

In the past decade it has become clear that super- 
massive black holes are a fundam ental ingredient of the 
universe (|Ferrarese fc Fordl 120051 ). Accretion onto their 
deep gravitational potential is responsible not only for 
some of the most powerful sources of light ([Lvnden-Belll 
1969), but also appears to be a key ingredient for the for- 
mation and evolut ion of galaxies (|Granato et al.l [2005 
iCroton et al.l [20061. Energy released by the accretion 
mechanism is believed to play a role in regulating the 
heating and cooling of interstellar gas and therefore the 
formation of stars. The "smoking gun" of this connection 
between galaxies and black holes is the tight correlation 
between the mass of black holes and the stellar velocity 
disper sion a* of their host galaxies, observed at small red - 
shifts (jFerrarese fc MerrittJ 120001: iGebhardt et al.l 120001) . 
This correlation represents the endpoint of the so-called 
coevolution of galaxies and black holes, though it is still 
unclear how galaxies and black holes coevolve across cos- 
mic time. Both galaxies and black holes are believed to 
grow by mergers and accretion from initial perturbations 
of the density field in the early universe; however, it is 
not yet known whether galaxies and black holes grow in 
lockstep, or one of the two forms first and subsequently 
acts as a seed for the other. 

The main challenge in mapping the coevolution across 
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cosmic time is determining black hole masses. Tradi- 
tional methods rely on spatially resolved kinematics of 
stars and gas within the gravitational sphere of influ- 
ence of the black hole. Therefore, they are only ap- 
plicab le with current technol ogy in the very local uni- 
verse (jFerrarese fc Ford 2005). Alternative methods are 
needed to measure black hole masses out to distances of 
several billion light years — lookback times correspond- 
ing to a sizable fraction of th e present-day age of the 
universe (jKomatsu et al.ll201lD . 

Reverberation mapping is the most promising method 
for measuring the masses of black holes powering ac- 
tive galactic nuclei (AGNs) at cosmologically interest- 
ing distances (f Blandfor d fc McKee]ll982t ; (Petersorill99l 
iPeterson et all 20041 ). The technique is made possible by 
the temporal variations in the intrinsic brightness of the 
central continuum source and by the subsequent response 
of line-emitting gas well within the gravitational sphere 
of influence of the black hole, known as the broad-line 
region (BLR). By measuring the time delay (or lag) r 
between the variations of the continuum and the varia- 
tions of the broad emission lines, the physical size of the 
BLR can be determined. In addition, the typical orbital 
velocity of the broad-line gas can be measured from the 
width of the broad lines themselves, 07. Combining this 
velocity measurement with the radius yields an estimate 
of the mass of the ce ntral black hole, Mbh= f&fcr/G 
([Peterson et al.ll2004[ ). where / is the so-called virial co- 
efficient, c is the speed of light, and G is the gravitational 
constant. 

Despite the simplicity of the reverberation mapping 
idea, its practical im plementation is beset with numer- 
ous difficulties (e.g.. iKrohkl 120011 ) . First, in its stan- 
dard implementation, the formula connecting black hole 
mass to spectral line width and the time lag includes a 
virial coefficient that depe nds on the unknown geometry 
of th e orbiting BLR gas (lOnken et all 1 20041: 1 Woo et all 
120101 : iGreene et afl 120101 iDecarli et all 120101 ) . Second, 
the adopted approach is usually indirect: the data are 



2 



Brewer et al. 



used to measure properties of the transfer function de- 
scribing the distribution of lags, or simply the mean lag. 
To constrain properties of the BLR gas distribution itself, 
another layer of modeling is necessary to reveal which 
possible BLR geo metries are consist ent with the inferred 
transfer function (jBentz et al.ll2010f ). 

We overcome these difficulties by applying a new 
framework to modeling high-quality spectrophotometric 
monitoring data of the Type 1 active nucleus of Arp 151. 
The combination of dynamical models with high-quality 
data allows us to characterize the structure of the BLR 
and achieve a direct determination of the mass of a su- 
permassive black hole using reverberation mapping. Our 
measurement is independent of any external information 
on the virial coefficient. 

2. DATA AND MODELING FRAMEWORK 

The data were collected as part of t he Lick AGN Mon- 
itoring Project (L AMP) campaign (jBentz et all 120091: 
IWalsh et al.1l2009l) . They consist of 84 epochs of photo- 
metric monitoring and 43 epochs of spectroscopic mon- 
itoring of the reg i on containing the H/3 emission line 
(jBentz et al.ll200a IWalsh et al.1 120091 ) . The continuum 
and line flux temporal series are shown in Figure [T] The 
H/3 spectral time series is illustrated in the left panel of 
Figure [5] (wavelength on the abscissa, epoch on the or- 
dinate; the epochs are approximately one day apart, but 
not exactly, due to gaps and scheduling issues). 

In order to infer the BLR geometry and the black hole 
mass we take a Bayesian In ference approach to the prob- 
lem Js ivia fc Skilling 200(f) . We construct a model char- 
acterized by a finite number of parameters describing the 
black hole mass, the spatial density profile of the BLR 
and its kinematic structure, and the intrinsic continuum 
light curve, denoting these parameters collectively by 3>. 
We define broad prior probability distributions for $ and 
then define the probability distribution for the data, con- 
ditional on knowledge of all of these properties p(D\<&). 
Given specific data, the prior distribution gets updated 
to the posterior distribution which describes knowledge 
of the parameters after taking into account the data, us- 
ing Bayes' rule p($\D = D*) tx p($)p(D\$)\ D = D , . In 
practice, we quantify our results by generating randomly 
sampled models from the posterior distributi on using a 
Nested Sampling method (Br ewer et al.|[20Toh ■ 

The physical model consists of a large number (1,000) 
of BLR clouds that are in orbit in the Keplerian poten- 
tial of the central black hole. The spatial distribution 
of the clouds is generated using a flexible model that is 
capable of representing generic geometries including thin 
disks and tori as well as complete spheres and shells. 
By applying spatially varying illumination to the clouds, 
we can also describe non-axisymmetric geometries. The 
cloud emission is assumed to respond linearly to con- 
tinuum variations. Therefore, the observed spectrum at 
a given time is the result of the continuum emission at 
earlier times, with time lags corresponding to the optical 
path from the central continuum to the cloud to the ob- 
server. The configuration of a model that well represents 
the data is illustrated in Figure [3l 

Fitting of models to the data requires knowledge of 
the continuum flux at all times, not just the measured 
times. To solve this problem, our method uses Gaussian 
processes to interpolate and extrapolate the continuum 
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Fig. 1. — Continuum flux time series from LAMP, and the corre- 
sponding H/3 flux time series. The curve drawn through the contin- 
uum data shows a realization of an interpolation of the continuum 
data using Gaussian processes. 

light curve taking errors into account. A typical intrinsic 
light curve generated by this process is shown in Fig- 
ure [TJ Thus, our results include uncertainty caused by 
the fact that we have noisy measurements of the contin- 
uum flux at a finite number of times. Our modeling of the 
continuum li ght curve is sim ilar to that independently 
developed bv lZu et"aL1 (|2010T ) and used outside of rever- 
beration mapping in studies of quasar variability (e.g., 
Kozlowski et all 120101 : Ikelly et al J 120091: iMacLeod et alJ 
2010ft . Add itional information on the method in general 
is given bv lPancoast et al.1 (|201lL hereafter Pll). 

3. THE DYNAMICAL MODEL 

We do not aim to infer the position and velocity of 
each cloud from the data, but rather to use the clouds 
to map out a spatial and dynamical model described by 
a small number of hyperparameters. This approach is 
equivalent to that followed by Pll, except that we are 
using Monte Carlo samples of clouds to represent the 
distribution of BLR gas as opposed to computing the 
density on a spatial grid. To generate a 3D distribution 
of BLR clouds, we start by generating an axisymmetric 
distribution in the x-y plane, and then apply rotations 
to "puff up" the model into a 3D configuration. Finally, 
we weight the clouds by a non-axisymmetric illumination 
function to model non-axisymmetric distributions of gas. 

The distance of a cloud from the black hole is pre- 
scribed according to r — F/i + (1 — F)Q, where F £ [0, 1] 
and Q is drawn from a gamma distribution with mean \i 
and standard deviation /3fi. With this prescription /i is 
the overall mean radius, F is the fraction of the mean 
radius that is due to the hard lower limit, and j3 6 [0, 1] 
describes the shape of the distribution; f3 w 1 is an ex- 
ponential distribution, and /3 « is a narrow normal 
distribution. The polar coordinate ^ of a cloud is chosen 
uniformly from [0, 2ir]. 

For the inference, the priors on these parameters are as 
follows. We use a uniform prior for /3, a scale-invariant cx 
1/x prior for /i (between generous limits) , and a uniform 
prior for black hole mass given fj,, such that the predicted 
line widths are on the order of those in the data, reducing 
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Fig. 2. — Top left: Measured spectrum of the broad H/3 emission 
line as a function of epoch; these are the data used for our inference. 
Top right: Model-predicted spectrum as a function of time, using 
parameter values chosen at random from the posterior distribution. 
The major features of the data (time variation, line widths) are 
reproduced by our model. Bottom: Example of the spectral line 
shape at two times, along with model fits. 



the volume of parameter space that needs to be explored. 

We then assign cloud velocities in a probabilistic man- 
ner (note that we can assign multiple velocities to each 
cloud in order to improve sampling of the phase space in 
a computationally efficient way; throughout this paper 
we adopt 100 velocities for each of the 1,000 clouds). As 
in PI 1, we assume that the only force acting on the BLR 
clouds is gravity from the central black hole. 

The total energy of a cloud at a distance r from the 
black hole, moving with angular momentum L, is given 

by 
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which has the minimum possible value 

m in 
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If we knew the position, energy, and angular momentum 
of a cloud, we could solve for the radial velocity, 
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We choose the negative (inbound) solution with proba- 
bility q and the outbound with probability 1 — q, a free 
parameter. For solutions to exist, the angular momen- 
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Circular orbits are obtained if we set the energy and 
angular momentum to 
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To get elliptical orbits, instead of assigning E and L the 
exact circular values above, we assign them at random 
from the following probability distributions: 
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l + exp(-x) / 
where x ~ A/"(0, A 2 ), 
and 



\L\ < L„ 
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These probability distributions are centered around the 
values for circular velocities, but the parameter A de- 
scribes the dispersion, or how noncircular typical orbits 
will be. The circular orbit formulae are reproduced when 
A — > 0. The probability distributions for E and L given 
three different values for A are shown in Figure HJ along 
with the corresponding line shapes. For the inference, 
we use a uniform prior on A between and 1. 

The main advantage of our implementation of this 
method with respect to that of Pll is that we can gen- 
erate relatively broad distributions of E and L with a 
smaller number of parameters. The drawback is that this 
model is technically nonstationary (in the case of noncir- 
cular orbits) and would change if it were allowed to evolve 
in a dynamically self-consistent way. Our model can be 
thought of as describing the time-invariant illuminated 
part of the full phase-space distribution, even though the 
underlying particles are actually flowing through the re- 
gion. Wc checked that this simplifying assumption does 
not bias our inference on the black hole mass by analysis 
of the data with the code developed by Pll. The results 
are consistent with the ones presented here. 

We then rotate the models (positions and velocities 
of the clouds) by an appropriate distribution of angles 
to generate an axisymmetric distribution of angular mo- 
mentum vectors. The first rotation is about the y axis 
by a small random angle; the typical size of these angles 
determines the opening angle of the disk or torus. We 
then rotate around the z axis by random angles to re- 
store the axisymmetry of the model. Finally, we rotate 
again about the y axis, by the inclination angle to model 
the inclination of the system with respect to the line of 
sight. For the inference, we use uniform priors on the 
opening angle and the inclination angle. 

To obtain non-axisymmetric models, in order to repro- 
duce the line asymmetry, we weight each cloud by a sim- 
ple spatially varying illumination function. In spherical 
polar coordinates, this function is 



1 



W(r, <f>,6) = - + k cos 4>, 



(11) 
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Fig. 3. — The distribution of extra path lengths the light must travel from the central engine to a BLR cloud and then to the observer is 
the cause of the delayed response of the emission-line flux, and the variations in line shape. The distribution of BLR gas in this diagram 
corresponds to a probable configuration inferred from the Arp 151 data. 




where k G [— |] is a free parameter with a uniform 
prior. Positive k illuminates the front portion of the 
BLR, negative k illuminates the back. We also imple- 
mented, as a secondary check, a model with a different 
(linear in x) functional form for the illumination, but 
this model did not reproduce the data as well, although 
the final black hole mass estimate was similar. In pre- 
dicting the observed spectra, we included the narrow-line 
component as a constant which does not respond to the 
continuum variations. 

4. RESULTS AND CONCLUSIONS 

Our results are presented in Figure El We find that the 
geometry of the BLR is well described by a thick disk 
or torus (opening angle 68.91^2 degrees), viewed at an 
inclination 67.8 ± 7.8 degrees (where 0° = edge-on and 
90° = face-on), as depicted in Figure |3l This geometry 
is consistent with the hypothesis that Type 1 AGNs are 
viewed close to face-on, if the dusty torus is coplanar with 
the BLR. The mean radius /i of the disk is 3.99 ± 1.25 
light days, and the radial profile is inferred to be close to 
exponential (j3 — O.SG^ig)- We note that four light days 
corresponds to < 10~ 5 arcsec at the distance of Arp 151 
(redshift 0.021), an angular size more than a thousand 
times smaller than what can be resolved even with the 
Hubble Space Telescope. The orbits of the BLR clouds 
are found to depart significantly from circular, and there 



is no strong evidence favoring either inflow or outflow 
(see Figure [5]). 

We note that the posterior distribution does not rule 
out radial profiles that peak close to r — 0, which is 
unphysic al due to the high ioniz ing flux from the accre- 
tion disk ([Korista fc Go ad 2004). However, conditioning 
on the peak of the radial profile being at the high end 
of the posterior does not significantly change any infer- 
ences except for the inflow fraction q: inflow becomes fa- 
vored by a ratio of 70:30 if we assume the density peaks 
at r > 1 light d ay. Thus, there is w eak evidence for in- 
flow, as found bv iBentz et al.l ()2010[ ). If inflow is present, 
then the front of the disk must be more visible than the 
back (i.e., k > 0), either because of obscuring material 
or nonuniform illumination effects. 

By marginalizing over all of the model parameters we 
derive the posterior probability distribution function for 
the central black hole mass. The median and 68% cred- 
ible interval are io 6 - 51±0 - 28 M©. This is lower than, but 
overlaps with, the value of io 6 - 85±0 - 07 M Q obtained by 
iBentz et al.l (2009) assuming log 10 / = 0.74 based on re- 
quiring active and inactive galaxies to obey the same 
correlation betwe en Mrh and host- galaxy stellar veloc- 
ity dispersion cr* (jOnken et al. 2004), and neglecting un- 
certainty in /. Recent measurements suggest that the 
intrinsic uncertainty in / from this method is at least 
0.4 dex (|Woo et al.'(l20la iGreene et aLll2010h . Revers- 
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Joint and marginal posterior distributions for the parameters of the BLR geometry and kinematics. The strongest correlation 
the inclination angle and the opening angle of the disk. Both of these parameters are also strongly correlated with the black 



ing the traditional argument, our measurement implies 
that log 10 / = 0.40 ±0.28, a low value, for this particular 
system. Modeling a larger sample of systems and com- 
paring the results of traditional methods with our direct 
approach would allow us to test the assumption that ac- 
tive galaxies obey the same Mbh-j* relation as inac- 
tive galaxies (iDavies et all 120061: Ifficks fc Malkanl 120081; 
lOnken et alj|2007j) . 

To summarize, our measurement has three key advan- 
tages with respect to traditional methods. First, it is 
direct, independent of any assumption regarding the cor- 
relations between supermassive black holes and their host 
galaxies, thus allowing us to test this assumption. Sec- 
ond, it provides a more precise measurement (smaller for- 
mal uncertainties) than traditional methods. Finally, the 
outputs of the inference procedure are physical properties 
of the BLR, rather than a transfer function, bypassing 
the need for an additional modeling step. 

We conclude by listing some of the limitations of this 
work and the prospects for future improvements. In our 
dynamical model we neglect radiation-pressure support 
and the dynamical influence of the accretion disk itself, 
as is commonly the case in traditional reverberation map- 
ping analysis. We also neglect the optical d epth of the 
BLR clouds themselves (jBottorff et al.lll997D . If the mo- 
tion of the BLR clouds were partially supported by ra- 
diation pressure, then we would be underestimating the 
mass of the central black hole. This is currently a topic 



of debate (|Marconi et alj2008t iNetzer fc Marzianil201(l , 
and external information needs to be used to break the 
degeneracy between black hole mass and pressure sup- 
port. Once external information is available, it can easily 
be taken into account in interpreting our results. 

On a more detailed level, our model for the spatial 
profile of the BLR is significantly oversimplified with re- 
spect to the real physical picture. Thus, we cannot all 
features of the observed line profiles (Figure [2|). Neglect- 
ing this systematic uncertainty would lead us to under- 
estimate the uncertainty in the parameter values. In this 
study, we addressed this issue by inflating the observed 
error bars (see Figure O until the model reproduced just 
the macroscopic features of the emission lines. In other 
words, we do not expect to be able to model all features 
of the lines, to within the given noise level. Inflating the 
measurement error bars on the data protects our results 
from some (but not all) systematic errors, particularly 
those that would result in fluctuations in t he line profiles 
smalle r than the domain of the data. See I Brewer et al.l 
(|201 If ) for a discussion of this issue. 

For these reasons, our uncertainty in the black hole 
mass is significantly larger than what the method can 
in principle deliver for data of comparable quality in the 
absence of modeling errors, ~ 0.05 dex (Pll). In ad- 
dition, our model neglects collisional effects as well as 
anisotropic winds, which could change somewhat the dy- 
namics and geometry of the BLR, but should not affect 
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the inference on black hole mass as long as gravity is the 
dominant force. 

Thus, modeling uncertainties dominate over measure- 
ment uncertainties, driving the total uncertainty in black 
hole mass. Therefore, the next step toward improving the 
overall precision of the measurement is to develop more 
flexible and physically realistic models. Such models will 
also allow us to explore in more detail the kinematics of 
the BLR. 
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